Finite size behaviors of critical Ising model on a rectangle with free boundaries 
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Using the bond-propagation algorithm, we study the Ising model on a rectangle of size M x N 
with free boundaries. For hve aspect ratios p = M/N — 1, 2, 4, 8, 16, the critical free energy, internal 
energy and specific heat are calculated. The largest size reached is M x N = 64 x 10 6 . The accuracy 
of the free energy reaches 10 -26 . Basing on these accurate data, we determine exact expansions of 
the critical free energy, internal energy and specific heat. With these expansions, we extract the 
bulk, surface and corner parts of free energy, internal energy and specific heat. The fitted bulk 
free energy density is given by /oo = 0.92969539834161021499(1), comparing with Onsager's exact 
result f x = 0.92969539834161021506 ■ ■ ■ . We prove the conformal field theory(CFT) prediction of 
the corner free energy, in which the central charge of the Ising model is found to be c = 0.5±1 x 10~ 10 
comparing with the CFT result c = 0.5. We find that not only the corner free energy but also the 
corner internal energy and specific heat are geometry independent, i.e., independent of aspect ratio. 
The implication of this finding on the finite scaling is discussed. In the second order correction 
of the free energy, we prove the geometry dependence predicted by CFT and find out a geometry 
independent constant beyond CFT. High order corrections are also obtained. 

PACS numbers: 75.10.Nr,02.70.-c, 05.50.+q, 75.10.Hk 



Introduction. Finite size effect has been attracting 
tremendous interest in the study of condensed matter 
physics. It becomes of practical interest due to the re- 
cent progresses in fine processing technologies, which has 
enabled the fabrication of nanoscale materials with novel 
shapes [H-Q . Exact solutions have been playing a key role 
in determining the form of finite size scaling. Ferdinand 
and Fisher Q pioneered on the 2D Ising model on a fi- 
nite size lattice, which extended Onsager's exact solution 
Q and stimulated the ideas of finite size scaling. Since 
then, exact results of the model on finite size lattices with 
various boundaries have been studied intensively (3- 12 1 . 
Detailed knowledge has been obtained for the torus case 
d, Q , for helical boundary condition , for Brascamp- 
Kunz boundary condition [ll, 12 1 and for infinitely long 
cylinder [l3| ■ However in the jigsaw puzzle of the 2D Ising 
model an important piece is still missed, which is the ex- 
act result on a rectangle with free boundaries. Generally 
speaking the rectangle geometry system is very interest- 
ing in its own right. It is, for instance, the natural one 
to consider in the case of quenches for one dimensional 
quantum systems with open boundaries. In 2D it is the 
simplest geometry to study transport properties for An- 
derson localization. Although there are Monte Carlo and 
transfer matrix studies on this problem 14, lij], the ac- 



curacy or the system sizes of the results are not enough 
to extract the finite size corrections. Meanwhile, for 2D 
critical systems, a huge amount of knowledge has been 
obtained by the application of the powerful techniques 
of integrability and conformal field theory (CFT) [16r - 
lilj ] . Cardy and Peschel predicted that the next subdom- 
inant contribution to the free energy on a square comes 
from the corners 17|, which is universal, and related to 
the central charge c in the continuum limit. Kleban and 



Vassileva [18j extended the study of the free energy on 
a rectangle. They further derived a geometry dependent 
term to the free energy. However, they did not determine 
a geometry independent additive constant in the coeffi- 
cient. Till now there is few evidence for these predictions 
from exact solutions or numerical calculations. Further- 
more, as far as we know, there is no detailed study of 
the internal energy and the specific heat on a rectangle 
neither by CFT nor by exact solutions/numerical calcu- 
lations. 

Recently an efficient bond propagation (BP) algorithm 
was developed for computing the partition function of 
the Ising model in two dimensions, which is exact to 
machine precision and works for any planar network of 
Ising spins with arbitrary bond strengths fljl . I20I ] . It is 
also much faster than Monte Carlo simulation, and costs 
quite moderate memory comparing with the transfer ma- 
trix method. Very large system size can be reached. The 
BP algorithm is thus a powerful tool to study the Ising 
model on a rectangle with free edges and corners. In 
this letter we apply the BP algorithm to study the Ising 
model on an M x N rectangle. We obtain finite size data 
of the critical free energy /, internal energy U and spe- 
cific heat c. By fitting these data we find that the exact 
expansion of the critical free energy, internal energy and 
specific heat can be written in the following form 
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where S = M x N is the area of the system. /oo,f7oo 
are the bulk term, / sur /, U sur f the surface coefficient, 
fcom,U corn the corner coefficient for the free energy, 
internal energy respectively. c sur f , c cor „ are the corre- 
sponding coefficients for the specific heat. We find the 
fitted values foo , f SU rf , Uoo , U sur f , An are excellently con- 
sistent to the exactly known results (#j8l. Hol . [ll| . The cor- 



ner free energy f corn — c/8 [ljlll8| is proved. More over 



the geometry independent constant, which is ignored in 
Kleban and Vassileva work [18|, is determined. We also 



find the corner contributions U corn and c c 



which are 



independent of aspect ratio. As far as we know, no previ- 
ous studies predict such terms. We start by introducing 
the BP algorithm briefly Then we present our numerical 
results and analysis. 

Method. The schematic of the BP algorithm is shown 
in Fig. [TJ In this algorithm BP series reduction, BP Y- 
A transformation and its inverse are the building blocks. 
By successively integrating in and then integrating out 
spin degrees of freedom in a way that only introduces lo- 
cal changes to the network, this algorithm progressively 
moves degrees of freedom to an open edge of the network, 
where they are eliminated. The transformation in each 
step is exact. The numerical accuracy is limited by ma- 
chine's precision, which is the round-off error 10~ 32 in 
the quadruple precision. The BP algorithm needs about 
N 3, steps to calculate the free energy of an N x N lat- 
tice (much faster than other numerical method). There- 
fore the total error is approximately iV 3 / 2 x 10~ 32 . This 
estimation has been verified in the following way: We 
compared the results obtained using double precision, in 
which there are 16 effective decimal digits, and those us- 
ing quadruple precision. Because the latter results are 
much more accurate than the formal, we can estimate the 
error in double precision results by taking the quadruple 
results as the exact results. We thus found that the error 
is about TV 3 / 2 x 10~ 16 . In our calculation, the largest 
size reached is M = N = 8000, the round-off error is less 
than 10 -26 . The partition function of the Ising model on 
a 2D square lattice is 



Z = ^2 exp(/3 ^2 
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where the nearest neighbor couplings are dimensionless 
and /3 is the inverse temperature. The free energy den- 
sity, internal energy per spin and specific heat density are 
defined by 
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(a) BP series 



(b) BP A-Y 
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FIG. 1. (Color online) (a) and (b) are building blocks of 
the BP algorithm: BP series and BP A — Y transformation, 
respectively, (c) The schematic of BP algorithm, (cl) The 
BP series is applied to the three spins and two bonds in red. 
(c2) The A — Y transformation is applied to the three spins 
and three bonds in red. (c6) Finally only two spins and one 
bond are left. 



respectively. With the BP algorithm, we get the free 
energy density / at the critical /3 C = \ ln(f + y/2) = 
0.44068679 • • ■ directly. The internal energy and specific 
heat are calculated by using a differentiation method 
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In our calculation, A/3 = 10 7 is used. The error in the 
calculation of U due to the finite A/3 is approximately 
i(A^) 2 9 3 //9/3 3 , which is less than fO" 11 for the lattice 
with N < 2000. The error in the calculation of c due 
to finite A/3 is i (A/3) 2 <9 4 //<9/3 4 , which is less than f0~ 9 
for N < 2000. Another source of error is the accumu- 
lated error due to round-off, which is less than 10 -26 . 
This error is amplified by l/A/3(l/A/3 2 ) times in the cal- 
culations of U (c), which is thus around fO -19 (10 -12 ). 
Therefore the errors of U, c are mainly caused by the finite 
A/3 in the differentiation. In other words the accuracies 
of the free energy, internal energy and specific heat are 
10~ 26 , KT U , 10~ 9 respectively. 

The calculations have been carried out for various as- 
pect ratios p = M/N = 1,2,4,8, 16 on an M x N rect- 
angular lattice. For M — N, the calculation was carried 
out from N = 30 to N = 8000. For p = 2,4,8,16, the 
calculated lattice sizes are from N = 30 to N — 2000. 

Critical free energy. We fit the data of free energy 
density with the formula given by Eq. (JT|) with k from 
1 to 8. The fitting method is the Levenberg-Marquardt 
method (2ll | for nonlinear fit. The standard deviation 

(SD) is defined by SD = y/^ifi - f^ fit) ) 2 /(n d - n f ) 
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with fi the numerical data, *^ the value given by the 
fitting formula, rid the data number and nj the number 
of fitting parameters. For all p, SD reaches 10 -20 . 

The high accuracy can be seen from the bulk 
value foo. For p = 1, the fitted values of is 
0.92969539834161021499(1). The asymptotic bulk value 
of the free energy density should be the same as the ex- 
act result given by OnsagerQ, i.e., foo = ln\/2 + \G = 
0.92969539834161021506 • • • , where G = 1-^ + ^-^ + 
• • • . Our estimation is consistent to it in the accuracy of 
10~ 20 . For other p the consistency is the same. 

According to finite size scaling, the surface correction 
term stems from free edges. This correction for the model 
on an infinitely long strip with two free edges has been 
found through exact solution to be (D\ — \ ln(l+Vz))/iV, 
where D x = /* ln[l + ^2(1 - cos6>) 1/2 (3 - cosfl)" 1 / 2 ] « 
0.2589553765253 % Note the infinitely long strip can 
be considered as a rectangle with a aspect ratio p — > oo. 
Thus the surface correction for a rectangle should be 
f surf (M + N)/S with f surf = D 1 - ±ln(l + V2) = 
-0.1817314169844187569 Our fitting results co- 

incide with this prediction precisely. For p — 1,2, 
the fitted values -f sur f are 0.18173141698441877(3), 
0.1817314169844188(1) respectively, and for p = 4, 8, 16, 
they are all 0.1817314169844188(2). The correction term 
fcom In S/S stems from the corners of the rectangle. Ac- 
cording to Cardy and Peschel [13], the four corners, at 
M = N, give rise to the term f \n(N)/N 2 , where c = 1/2 
is the central charge. Kleban and Vassileva [l~8| extended 
the results to a rectangle with M 7^ N using CFT: The 
corner free energy is chi(S)/(8S), which yields 



fc 



= -= 0.0625 



(7) 



In our fitting f corn are found to be 0.0625+ (2±3) x 10 14 , 
0.0625 + (1 ± 2) x 10" 12 ,0.0625 + (1 ± 5) x 10- 12 ,0.0625 + 
(1 ± 5) x 10~ 12 ,0.0625 ± 1 x 10" 11 for p = 1, 2, 4, 8, 16 re- 
spectively. All of them lead to the central charge c = 0.5 
with the error less than 10~ 10 . There is a 1/iV 2 correc- 
tion for the infinitely long strip with the coefficient 
ctt/24 « 0.0654498. Kleban and Vassileva [H proved 
that, for a finite rectangle, the correction is proportional 
to 1/5 = 1/(M x N) with the coefficient 



(8) 



where rj(q) = q 1 / 24 



-2-n-p 



= e ""r^q- = 
ctt/24: which 



Z 1/24 nXl(l-0 with q 
e -27r/p. At the limit p ->• 00, A\p~ x - 
recovers the infinitely long strip result. 

However, Kleban and Vassileva mentioned that, in 
their derivation, a possible geometry-independent addi- 
tive constant was ignored [18j. In other words, the coef- 
ficient of 1/5 should be 



where Fq is the constant which contributes Fq/S — 
Fo/(pN 2 ) to the free energy density, which tends to 
in the infinitely long strip limit. For p = 1,2,4,8,16, 
the values of A\ are 0.065918017562, 0.087578866954, 
0.175155990232, 0.393633679243, 0.873910756056 re- 
spectively. By comparing the fitted A\ (see Tab. Q} and 
the theoretical A[ , we obtained this constant 



F = -0.0049488147(2) 



(10) 



The other parameters A2 , • • • , As are fitted and listed 
in Tab. |H For the infinitely long strip, A 2 p~ 3/2 , A 3 p~ 2 
and Aip~ b / 2 , at the limit p 00, correspond to the 
coefficients of N~ 3 ,N~ 4 and N~ 5 terms in the finite 
size free energy expansion respectively, which have been 
obtained as -0.04616(2), 0.024(1), 0.69(6), by using nu- 
merical transfer matrix techniques [22l |. Simple extrap- 
olations of the fitted values of A2 , A3 , A4 show that the 
magnitude of A2 , A3 agree with the transfer matrix re- 
sults, but A4 does not. We have tried other forms of 
formula to fit the critical free energy data. For example, 
we added the terms ln5/5 3 / 2 ,ln5/5 2 in the fitting for- 
mula and found that the corresponding coefficients are 
extremely small (less than 10 -7 ). We concluded that the 
logarithmic correction only appears in the corner term 
In S/S. We note here that, in the asymptotic expan- 
sion of the free energy for the Ising model with periodic 
boundary conditions [8|] , with Braskamp-Kunz boun dary 
conditions [ll[ and with helical boundary conditions [ldj . 
only integer powers of S appear. 

Critical internal energy. We fit the data of critical in- 
ternal energy with the formula given by Eq. Q with 
k from 1 to 4. The bulk value Uoo is known to be 
V2 « 1.41421356237 @. Our fit of is 1.4142135624 
for the five aspect ratios. The leading correction 

should be caused by the edges. In the exact result of 
Au-Yang and Fisher 0] on the strip with two free edges, 
the edges' correction is given by — ^ hi N/N. We con- 
jecture that, on the rectangle with four free edges, the 
edge correction is given by - 1 (M In N + N In M)/(MN) 
with the coefficient U sur f(p) = — § = —0.636619773 
This conjecture is proved by the fitted U sur f for all p 
(U sur f = —0.6366198...), which agree with the predicted 
value in the accuracy 10 -7 . 

The term B\ / S 1 / 2 is in fact scaled as 1 /N . In the infin- 
ity long strip limit, the coefficient before 1/N is known 
as |(|ln2 + 7 - f - hiTr) « 0.683158 0. As we can 
see in Tab. UH B\j ^fp indeed approaches this limit as p 
increases. 

Following the convention in the critical free energy, we 
write the coefficient of (In S)/S as U corn . We found 



U corn {p) w -0.45025(3), 



(11) 



A x = A[+F , 



(9) 



which is independent of aspect ratio p. Apparently, this 
contribution becomes zero in the limit of infinitely long 
strip, in which there is no corner. From this point, it is 
rational to call this term corner's correction. 
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TABLE I. The fitted parameters of Eq. |T]) for the critical free energy. 



parameter 


p=l 


p = 2 


p = 4 


p = 8 


p = 16 


Ai 


0.060969202833(3) 


0.08263005222(2) 


0.170207175502(4) 


0.388684864512(8) 


0.86896194132(2) 


A 2 


0.0883883476(2) 


0.0595556310(1) 


-0.1009034881(4) 


-0.666287857(1) 


-2.423248445(3) 


A 3 


-0.0175362651(2) 


0.04067435(2) 


0.41842401(7) 


2.1854954(3) 


9.76558(1) 


M 


-0.02405666(2) 


-0.174428(2) 


-1.290133(9) 


-8.15503(5) 


-48.556(3) 


A 5 


0.066893(8) 


0.4351(1) 


3.9657(8) 


33.676(6) 


277.21(5) 


A 6 


-0.16558(2) 


-1.147(4) 


-13.85(5) 


-161.7(5) 


-1857(6) 


A r 


0.3851(4) 


3.22(8) 


53(1) 


86(2)0 


137(4)00 


A s 


-0.703(3) 


-7.6(8) 


-17(2)0 


-387(4)0 


-8(1)0000 



TABLE II. The fitted parameters of Eq. ([2]) for the critical internal energy per spin. 



parameter p — 1 p — 2 p = 4 p= 8 p = 16 

U surf -0.63661981(1) -0.63661983(3) -0.63661985(4) -0.63661987(5) -0.63661988(5) 

Bi -0.1213621(2) -0.2586182(5) -0.6453965(8) -1.266482(1) -2.151535(2) 

U CO m -0.450170(3) -0.45018(1) -0.45022(2) -0.45028(4) -0.45042(9) 

B 2 -0.98604(2) -1.1278(9) -1.7138(2) -3.2016(4) -6.4893(9) 

B 3 -0.2969(2) -0.101(1) 0.833(3) 3.943(8) 13.40(2) 

B4, -0.040(2) -0.41(1) -2.40(4) -11.2(2) -47.7(7) 



The other parameters B\ , B2 , -B3 , B4 are fitted and 
listed in Tab. [TTJ Again we have tried other forms of 
formula to fit the critical internal energy. The terms 
lnSyS 13 / 2 , In5/S* 2 are excluded considering the coeffi- 
cients are extremely small. Moreover the standard de- 
viations of the fits with these terms are much larger than 
those without them. 

Critical specific heat. The data of the critical specific 
heat are fitted using the formula given by Eq. ([3]) with 
k from 1 to 4. The leading term Aq hiN is known from 
Onsager's exact result [Bj], which reads A = -[ln(l + 
v^)] 2 ~ 0.494538589. Our fitting gives A « 0.49453858. 
The other fitted parameters are shown in Tab. Mil The 
coefficient cq increases with the aspect ratio. For the 
strip case it is known that — co(p = 00) = (| In 2 + 7 — 
14£(3)/tt 2 - I -hiTr) w 0.3125538 0. c approaches this 
limit as p — > 00 obviously, see Tab. IIII1 However we have 
not obtained an analytical expression for the dependence 
of Co on p. 

The term (M In N + N In M)/S is the next order cor- 
rection. Its coefficient c sur f is independent of p, and its 
average value over p is c sur f — 0.524529(3). Note that 
this term is absent in the torus case and not men- 
tioned in the long strip case Q , but exists in the cylinder 
case with Brascamp-Kunz boundary conditions (ill [l2j ■ 

The corner term c corn also seems independent of aspect 
ratio p and its average value is about 



= 0.368(1). 



(12) 



the terms (lnS*) 3 /^, (lnS) 2 /S in the fitting formula and 
found that their coefficients are extremely small. 

Conclusion. Using the BP algorithm, we studied the 
Ising model on a rectangle of size MxN with free bound- 
aries. For five aspect ratios p — 1, 2, 4, 8, 16, the critical 
free energy, internal energy and specific heat were calcu- 
lated. We brought numerical evidence that the correction 
term f corn stems from the corners of the rectangle is in- 
deed universal and is proportional to the central charge 
c. We also found that the terms U corn and c corn are in- 
dependent from the aspect ratio p. In order to check 
whether or not the terms U corn and c corn are univer- 
sal quantities, it is useful to extend our study to the 
Ising model on other types of lattices, e.g., the triangu- 
lar, honeycomb, Kagome lattices. In such studies, the 
BP algorithm which is suitable for any planar network 
of Ising spins with arbitrary bond strengths [ill [2(| is 
still a powerful tool. In addition, we can enhance the 
accuracy of the internal energy and specific heat by us- 
ing an extended BP algorithm to calculate the internal 
energy without the using of differentiation 20] . Moreover 
the finite size effects of the correlation can be investi- 



We have also tried many other forms of formula to fit 
the critical specific heat data. For example, we added 



gated by using the site propagation algorithm [23j. As 
we have shown, the sharp corners induce remarkable ef- 
fects in critical region not only on the free energy, but 
also on the internal energy and the specific heat. It is 
expected that the sharp corners can induce remarkable 
effects on other properties of finite size systems in critical 
regime, for example, thermal conductivity, electric con- 
ductivity, etc. All of these effects should be observable 
in experiments. 



5 



TABLE III. The fitted parameters of Eq. ([3]) for the critical specific heat. 



parameter p — 1 p — 2 p = 4 P — 8 p = 16 

c -0.57078599(3) -0.44276294(2) -0.37766115(2) -0.34510742(2) -0.32883055(2) 

c surf 0.524516(2) 0.524533(2) 0.524531(2) 0.524530(3) 0.524529(3) 

D a -0.34928(3) -0.29524(1) -0.17826(1) -0.05986(2) 0.03117(2) 

c corn 0.3683(7) 0.3704(3) 0.3696(7) 0.368(1) 0.365(3) 

D 2 1.144(5) 1.293(2) 1.914(5) 3.42(1) 6.67(3) 

D 3 0.02(4) -0.13(2) -1.25(7) -4.8(2) -15.4(6) 

D 4 0.8(3) 0.7(2) 3.6(8) 16(3) 7(1)0 
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